Integrability of an extended d + id-wacve pairing Hamiltonian 
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We introduce an integrable Hamiltonian which is an extended d + id-wave pairing model. The 
integrability is deduced from a duality relation with the Richardson-Gaudin (s-wave) pairing model, 
and associated to this there exists an exact Bethe ansatz solution. We study this system using 
the continuum limit approach and solve the corresponding singular integral equation obtained from 
the Bethe ansatz solution. We also conduct a mean-field analysis and show that results from these 
two approaches coincide for the ground state in the continuum limit. We identify instances of the 
integrable system where the excitation spectrum is gapless, and discuss connections to non-integrable 
models with d + id-wave pairing interactions through the mean-field analysis. 
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I. INTRODUCTION 

The theory of Bardeen, Cooper and SchriefTer 1 (BCS) 
to describe the superconducting properties of metallic 
materials is now more than fifty years old. An exact 
Bethe ansatz solution in the case of a reduced BCS 
Hamiltonian with uniform coupling parameters was ob- 
tained by Richardson and Sherman-^. This same system 
was also studied by Gaudin using an electrostatic anal- 
ogy and the continuum limit approach. 4 Today it is of- 
ten referred to in the literature as the Richardson model, 
Richardson-Gaudin model, or s-wave model. Some time 
later it was shown that this system is integrable^ and 
following on from this the integrability and exact so- 
lution were unified through approaches involving the 
quantum inverse scattering method based on the Yang- 
Baxter equation,— and Gaudin algebra methods^ From 
those developments up to the present day, the s-wave 
pairing model has continued to be the subject of many 
investigations including the implementation of numeri- 
cal techniques for solving the associated Bethe ansatz 
equations^— studies using analytic methods^ - — the 
computation of correlation functions^ - — and quantum 
dynamics , 17 ' 18 and relations to conformal field theory.™ 

There have also been several searches for examples 
of integrable, exactly solvable BCS Hamiltonians which 
go beyond the s-wave case. Some examples were found 
whereby integrability was maintained by the inclusion 
of additional interaction terms besides the familiar pair- 
scattering interactioni22r— Recently, it was also found 
that integrability holds in the case of p + ip-wave 
pairing!^ - — An interesting feature of the p + ip-wave 
pairing case is that the model has a non-trivial ground- 
state phase diagram, in contrast to the s-wave model. A 
natural question to ask is whether the next in the se- 
quence, a Hamiltonian with d + id- wave pairing ; 27 ' 28 also 
admits an exact solution and if so what is the predicted 
phase diagram? Within the typical (Lie algebraic) ansatz 
for the model's eigenfunctions, a recent study 29 indicates 
that the answer is negative regarding the existence of an 
exact solution. Here we will show that with the inclusion 
of additional interactions, of the same quadratic number 



operator type as those encountered in earlier studies* 2 ^— 
an exact solution can be obtained. Our approach uses the 
same conserved operators as those of the s-wave model, 
but the Hamiltonian is constructed via a different combi- 
nation. We examine the continuum limit of the Bethe 
ansatz equations, and the gap and chemical potential 
equations from mean-field theory for the ground state, 
and find that these agree. The mean-field equations point 
to the existence of a critical point where the excitation 
spectrum is gapless. This critical point manifests itself 
in the Bethe ansatz equations as the point at which all 
ground-state roots are equal to zero. We finally discuss 
connections to non-integrable models with d + id-wave 
pairing interactions through the mean-field analysis. 

In Sect. 2 we establish integrability of the extended 
d + id- wave pairing Hamiltonian by using a duality trans- 
formation applied to the s-wave pairing Hamiltonian con- 
served operators. In Sect. 3 we use the continuum limit 
of the Bethe ansatz equations to study the the system by 
solving the corresponding singular integral equation. We 
also discuss the arc in the complex plane that contains 
the roots of the Bethe ansatz equations in the contin- 
uum limit for uniformly distributed single particle levels 
(free fermions in two dimensions). In Sect. 4, we use 
the mean-field approximation to obtain the chemical po- 
tential and gap equations. We show that these results 
coincide with results of Sect. 3 in the continuum limit, 
and conclude that the results found in Sect. 3 refer to 
the ground-state of the system. We also discuss connec- 
tions to critical points for non-integrable cases through 
the mean-field analysis. Concluding remarks are given in 
Sect. 4. 



II. THE HAMILTONIAN AND 
INTEGRABILITY 

We first introduce the Hamiltonian of the singlet- 
pairing BCS model with d + id-wave symmetry. We 
take the canonical (i.e. particle number preserving) 
Hamiltonian whose mean-field approximation leads to 
the Bogoliubov-de Gennes equations with order param- 
eter having d + id-wave symmetry (e.g. see Read and 
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Green^Z, Sato et al<2£) up to quadratic approximation. 
Letting c k(T , ct denote annihilation and creation oper- 
ators for two-dimensional fermions of mass m with mo- 
mentum k = (k x ,k y ), the Hamiltonian reads 



(1) 



- G^(fc x + ik y ) 2 (k' x - ikyf^ci^c-k^ck'- 



The s-wave pairing Hamiltonian is given in terms of 
hardcore boson operators by 



3 = 1 j,k=l 



(4) 



k,k' 



The integrability of this Hamiltonian was shown by ex- 
plicit construction 5 of a set of commutating operators 
(integrals of motion). These operators read 



where a I denote the spin labels and G is the cou- 
pling constant (positive for an attractive interaction). It 
is more convenient to introduce the following Cooper pair 
operators — c^c^^, 6k = c_ k ^c k -|-, and Nk — & k ^& k . 
For any unpaired fermionic state the action of the pairing 
interaction is zero and we can decouple the Hilbert space 
into a product of paired and unpaired fermions states, for 
which the action of the Hamiltonian on the space for the 
unpaired fermions is automatically diagonal in the nat- 
ural basis. We can therefore exclude unpaired fermions. 
We consider m = 1, Sk = |k| 2 and k x +ik y = |k| exp(i0 k ) 
in which case the Hamiltonian takes the form 



H = ^£k^Vk-G^£ k £k' exp(2i</) k )exp(-2j0 k /)& k t 6 k /. 

k k,k' 

We next perform the following unitary transformation 

bk f = exp(-2i(f>k)bl, 6 k = exp(2i0 k )6 k . 

Using the integers rather than k to enumerate the mo- 
mentum states, and adding a constant term, the Hamil- 
tonian takes the form 



G L 



L L 

H = J2 ekNj -GJ2 zM b » + 5 E E ) 

3 = 1 j,k=l 3 = 1 



x>^ - f e e ^ { b \ bk + wt 



G 

3 = 1 ~ j,k=l 



+ 26 ik NjN k 



where the hardcore boson operators satisfy the following 
commutation relations 

h, b\] = Mi - 2&&), hM = [b],b\] = o 

as well as (&]) 2 = 0, Nf = Nj. Setting G = G/2 we now 
introduce an extended d + id-wave pairing Hamiltonian 

L L 

H = E e i N i ~ G E £ Mb]b k + bjbt + 2NjN k ), (3) 

3=1 j,k 

We will establish integrability of ((3]) using the integrabil- 
ity properties of the s-wave model and a duality transfor- 
mation involving the coupling constant G and the single 
particle energy levels e k . We assume that the e k > are 
ordered such that Sj < e k whenever j < k, but otherwise 
there are no constraints imposed. 



^ = -j7(2^--0+E- 



*7* 



(5) 



with 



9 = b 1 ® b + b ® b ] + i(27V - I) ® (2N - I). 

The eigenvalues of these operators (excluding the space 
of unpaired fermions) are given by^ 



1 I 

G 



' 9 ^ u - e k ^ e,- - v k ' 

k^j 3 1=1 3 



(6) 



where the v k are solutions of the corresponding Bethe 
ansatz equations 



2 ^ 1 ™ 1 

k=l j^l J 



(7) 



The Hamiltonian given by Eq. can be written in terms 
of these commutating operators (up to a additive con- 
stant) as 



H = -G'J2 



(GO 



3=1 



E 

j,k=l 



TjT k 



y\2 L Qf L 



with C = tfb + btf + 2{2N - I) 2 . The energy spectrum 
of dU) is given by 



M 



E = 2^, 



3=1 

where M represents the total number of Cooper pairs. 
The above discussion outlines the essential features 
concerning the integrability and exact solution of the 
Richardson-Gaudin model ([4]). 

We now consider the following change of variables and 
insert this into the expressions for the integrals of motion, 
their eigenvalue formulae, and the Bethe ansatz equa- 
tions: 



e -> e \ v -> y 1 . 



(9) 
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We have from Eqs. d5j) , (JSj) , and ([7]) respectively 

L 

a 7v, - n -I- ' 

G' 



1 1 \ - 



e k - ej 

£j£k \ - yi£j 



k+j 

L M 



k+j 



£fc 



E 



VI - £ 3 



2 x 



(10) 

(11) 
(12) 



Any combination of the conserved operators will com- 
mute with the full set of conserved oeprators and can 
be taken as an integrable Hamiltonian. Instead of the 
combination given by Eq. © that corresponds to the 
s-wave Hamiltonian ((3]), we take the following combina- 
tion, written in terms of the new parameters Ei given by 
Eq. © 



j=i 



(13) 



Next we implement the following transformation of the 
coupling constant 



4G- 1 = -G'- 1 + 2j2?j- 

We are now in a position to express the Hamiltonian of 
the extended d+id-w&ve pairing Hamiltonian in terms 
of the Hamiltonian H given by (|T3|) : 



L L 



yE^ + oEE^-E 



3=1 



3=1 k+j 



3=1 



The energy of ([3]) is found to be given by the roots of the 
Bethe ansatz equations 



1 L F 2 

— + Y ii 

2Gy; j-[yi{yi-Ei) 



M 
3+1 



V3 



vi - y-j 



l = 1,...,M 



through the expression 

M M M 



(14) 



1=1 



, ViVl ~ /_^~y 

1=1 3+1 3=1 



Eqs. (fT4| and (fl5jl provide the exact solution of the ex- 
tended d + id-wave pairing model ([3]) on the space with 
no unpaired fermions. Extending this solution to the full 
Hilbert space of states is straightforward by omitting the 
appropriate blocked levels from the summations over Sj, 
and by accounting for the contribution of the blocked 
states to the total energy. 



The eigenstates of Hamiltonian are also obtained in 
the above method. In terms of the solutions of (IT41 they 
read 



A I 



l*} = Il C fe)|0), C{y)=Y J 



k=l 



y~£k 



where |0) is the vacuum state. It can also be shown via 
a direct algebraic calculation that the following identity 
holds: 



[H,C M (Q)) =MC M - 1 (0)Qt(2G^£i - 1). 



This property is reminiscent of the p + ip-wave pairing 
modelj^— Given any eigenstate |\&) with energy E, the 
dressed state \^f') = C M (0) \^) is also an eigenstate with 
energy E whenever 2G^ /=1 £; = 1. In the following 
analysis it will be seen that states which are dressings of 
the vacuum state will play a distinguished role. 



III. CONTINUUM LIMIT AND SINGULAR 
INTEGRAL EQUATION 

In this Section we study the continuum limit of the 
Bethe ansatz equations, and look to solve the associated 
singular integral equation. The Bethe ansatz equations 
as given by (fT4|) and can be written as 



1 



G D ^ 

— + -2 + > 

yi yf friw- ^ 



M 



= 2£ 



i 



3+1 



yi - Vi 



i = l,...,Af 



for G = 2(M-1)- Land 



D 



1 

2G 



E< 



The Bethe ansatz equations for D = would ad- 
mit an electrostatic analogy as has been used in other 
studies i^£ ' 12 i 13 ' 16 i 25 ' 26 In that setting the solutions yj are 
equilibrium positions of M mobile +1 charges subject 
to electrostatic forces from —1/2 charges associated to 
the L fixed positions and a charge — C/2 at the ori- 
gin. Even though the electrostatic analogy no longer 
holds for D ^ 0, the equation can still be solved by 
use of complex analysis techniques. In the continuum 
limit G~\ M, L ->• oo with x = M/L and g = GL held 
constant, we assume that the roots yj are distributed 
according to a density function r(y) with support on an 
arc r in the complex plane. Introducing a second density 
function p(e) with support on an interval O = (0, lj) of 
the real line such that 



de p{e) = 1 



4 



we can write 



c d 



ds 



2/ y Jn £ - y 

where c = 2x — 1 and 



d = — — / dssp(s). 
2.9 7n 



0, (16) 



(17) 



We first consider the case T = Ta U Tb where Ta — 
(0, £a)i £a < w and Tb is an arc, cutting the real axis 
at ea, which is invariant under complex conjugation and 
has end points a, b such that a = b* = e + iS. In this case 
Eq. ([TO]) becomes, Vy G T B , 



' " 



fie 



y y Jo £ -v J SA £ -v 

j Mv') 
'v., v-y r 



p I W 

Jr R 



(18) 



We introduce a function h(y) which is analytic outside 
of 0, and Tb, with a branch cut on Tb such that 
2mr(y)\dy\ = {h + (y) - h_(y))dy Vy G T B . Specifically 



h £ -y y y 



R{y) = ^{y-a){y-b) 



(19) 
(20) 



where R(s) is positive on (sa,w) and negative on (0,£a)- 
Letting Cb denote a closed contour around Tb, the 
Cauchy principal value of the expression in Eq. (|18[) can 
be written as 



P[ 2\dy'\^ 1 =<f 

Jtb y-y Jc 



dy' h(y>) 

c B 2ni y'-y' 



(21) 



Using residues to evaluate (l2~Tj) for the ansatz (|T9|) . and 
substituting into (fTS)) leads to 



Pie) 



c d f eA , p(e) y , 

y y jo e-y ./^ e-y 



(e)iZ(e) 



e-y Jo £ -y Jo 



ds- 



T 1-™ \+U =0. 



\_ y y \2Vaby y z 
The above equation is satisfied by choosing 



P(e) 
\R(s) 



u = j, _ c , d{a + b) 



and 



ab 2abVab 
c d(a + b) 



ds4>(s) - — — 
o V ab 2ab\f ab 



(22) 



Given this solution we may now calculate quantities of 
interest. The first example is the filling fraction x: 



\dy\r{y) + / \dy\r{y) 
r A Jt b 

= Jo d£P<y£ ^ + ^iJ c dyk ^ 

1 c d 1 f u , . 
= 2 + 2 + ^o-2Jo d££m - 

However, as the parameter c also depends on the filling 
fraction, the x terms cancel and we obtain the following 
equation 



- = 2 I dssp(s) + 2\/ab dss<t>{s) 
9 Jo Jo 



(23) 



Now combining (|17l22l23p gives the following equation 
for the filling fraction 



1 ^fab 



dscj>{s) - (—-3. I dse(f>(e). (24) 



Wab 



Defining the intensive energy e as 

e = lim — 

we have from (|T5|l that e = X — 2gl 2 where 

1=1 \dy\yr{y) + / \dy\yr(y) 
Jr A Jt b 



f I dy 



d(^-Vb) 2 



Wab 



ds ep(e) 



i r 

— — J dss(f>(s)(2s — a — b) 
1 f u 1 f u 



d££(j)(£) 



(25) 



where in the last step we use (|17I23|) . 

Repeating the above calculations for T = Tb cutting 
the real axis at sa < 0, again invariant under complex 
conjugation and having end points a, b such that a — 
b* = e + iS leads to exactly the same Eqs. (|23|24I25|) . 



A. Equation of the arc 

We now specialise to the case of energy levels s uni- 
formly distributed in the interval [0, 1] (i.e. taking uj — 
1), such that p(s) = 1. This case corresponds to the dis- 
tribution of two-dimensional free fermions. Integrating 



5 



Eqs. (|23I24|) leads to 



- = 1 + 2^e 2 + <$ 2 (V(e - I) 2 + S 2 - Ve 2 + 5 2 ) 



2e\/e 2 + S 2 In 



v/(l - e) 2 + <5 2 + 1 - 



V^TT 2 



1 e / /(e-l) 2 + 5 2 
• T= 2 + 2 e 2 + * 2 



5 2 



2Ve 2 + S 2 



In 



v / (e-l) 2 + <5 2 + l-e 



Ve 2 + <5 2 - e 



For a given value of the coupling constant g and filling 
fraction x we can numerically obtain the end points e±i5. 
From these end points we can numerically evaluate the 
following integral that gives the equation of Tb in the 
complex plane: 



Re 



dy'h(y'] 



0, yer B . 



(26) 



For the uniform distribution p(e) = 1 it is easily ver- 
ified that Eq. (|23|) is satisfied by the choice g = 1, e = 
6 = giving e = X = 0. We recognise the corresponding 
eigenstate as the dressing of the vaccum as discussed in 
Sect. |TTJ Fig. 1 shows the limiting behaviour of the arc 
r as g — > 1, confirming this conclusion. We highlight 
that this scenario is very different from the case of the 
p+zp-wave model. In that model roots collapse at the ori- 
gin when the coupling crosses the so-called Moore-Read 
line for any finite number of particles, but the limiting 
behaviour of the arc in the continuum limit is a closed 
curve, of finite length, encircling the origi m 25 ' 26 

In the above calculations we have obtained a unique 
solution of the Bethe ansatz equations in the continuum 
limit under the assumption that the roots become dense 
on a connected curve in the complex plane. Now we 
will argue that this solution corresponds to the ground 
state of the Hamiltonian, by making comparison between 
the exact solution and results obtained from a mean-field 
approximation of the Hamiltonian. Using a mean-field 
approximation will also allow us to make some compar- 
ison between the properties of the d + id-wave pairing 
Hamiltonian (JTJ) and those of the exactly solvable ex- 
tended d + id-wave pairing Hamiltonian ([3]). 



IV. MEAN-FIELD APPROXIMATION 

We consider a Hamiltonian with the compactly written 
form 

H = H - GQ^Q - GQQ f - 2FH 2 , (27) 

L L 

where H = Y^e&iVfc, Q = ^fffcfefc and F, G > are 

fc=l k=l 

the two coupling parameters. Using a mean-field approx- 
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FIG. 1: Arcs representing solutions of the Bethe ansatz equa- 
tions in the continuum limit. Results are shown in the case 
of p(e) = l,x= 1/2 for (a) the arc T with g = 203/200, g = 
202/200 and g = 201/200; (b) the arc F B with g = 199/200, 
g = 198/200 and g = 197/200. The limiting behaviour indi- 
cates that the arc T contracts to a point at the origin when 
g = 1. This indicates that the corresponding eigenstate at 
g = 1 is the state obtained by dressing the vacuum. 



imation the Hamiltonian given by Eq. (|27|) becomes 
U = H 



AFT H + 2FT 2 - fi(N - M), 



(28) 



where A = 4G (Q), T = (H ), A = |A|, N is the Cooper 
pair number operator, M = (N), and fj, is a Lagrange 
multiplier which is introduced since the mean-field ap- 
proximation does not conserve particle number. We will 
assume that the total fermion number is even, so all 
fermions are paired in the ground state. Setting F = 
gives a mean-field approximation for ([l} (up to a canon- 
ical transformation), while setting F = G gives a mean- 
field approximation for ([3]) . Keeping F as a free variable 



6 



allows us to interpolate between these two cases. 
Setting s = 1 - 4JT and 



£{e 3 ) = ^{sE j - l ,y+E^ (29) 
the ground-state energy is found to be 



E = \^{se ] -y)-\Y J £{e 3 ) + ^+2FT 2 +iiM (30) 

3=1 J =1 

associated to the mean-field ground state 

L 

i* m f)=n(v+^])io) 



3=1 



where 



1 



SBj — fl 



i 



SBj — fl 



The elementary excitation energies are given by £(Ej)/2. 
It is seen that the excitation spectrum becomes gapless 
in the limit E\ — > when /i = 0. 

Through use of the Hellmann-Feynman theorem we 
may take partial derivatives of (|28p and (|30|) to gener- 
ate the following constraint equations: 



(31) 
(32) 

(33) 





1 

2G 


L 

= E 

3 = 1 
L 


L - 


- 2M 


= E 

3=1 


L 




L 


E £ 3 

3=1 


- 2T 


=E 

3=1 



Combining Eqs. (|31l33p gives 
1 + 4(G - F) 



2G 



(34) 



Equivalently, the explicit expression for ,s in term of i* 1 
and r leads to 



L L 

1 — s s \ -s \ - e,- 



2F 2G 



E £ 3+mE 



3 = 1 



(35) 



Using (|31I32I33[) the energy expression (j3"0")) may be 

rewritten as 



£■ = T - 2FT 1 



A 2 

8G 



(36) 



which is consistent with (1271) . 

First we consider the integrable case F — G for which 
(l35l) reduces to 



L L 

-o=E^+/ i E 



2G 



3=1 



3=1 



(37) 



Expressing (l29l) as £(e) = Vs 2 + A 2 y / (er — a)(e — b) we 
have 



Identifying 



a = b* 



Sfl 



s + iA 
s 2 + A 2 



+ A 2 ' 



<5 = 



A/i 



A 2: 



in the continuum limit we find that (|37|32|) and e calcu- 
lated through (|3"6")) coincide with the exact solution ex- 
pressions (|23I24[) and e calculated through (12B1 respec- 
tively. This gives confidence that the mean-field approxi- 
mation is valid when F = G, while simultaneously show- 
ing that the solution found in the continuum limit in 
Sect. [ni]does correspond to the ground state of the sys- 
tem. Gapless excitations in the mean-field analysis arise 
when /i — which, from (|37l) . corresponds to 



1 

2G 



L 

3 = 1 



independent of x. The above relation is exactly the con- 
dition for dressing states, presented in Sect. [TTJ 

Finally we consider the mean-field equations for gen- 
eral F 7^ G when n — 0. Here we deduce from (|31I32I33[) 
that the excitation spectrum is gapless when 

L 

2G£>-1 

3 = 1 



4(G-F)]T £ , 

3 = 1 

Qualitatively, the d + id- wave case where F = and the 
extended d + id-wave case with F = G are very similar at 
the level of the mean-field analysis. For fixed x, in both 
instances there is at most one value of the coupling G 
at which the spectrum may be gapless. One significant 
difference however is that gapless excitations for the case 
F = can only occur when x < 1/2, whereas for F = G 
gapless excitations occur for all filling fractions. 
For non-zero £j we find when fi = 



In 



Vs 2 + A 2 



Vs 2 + A 5 



independent of j. Alternatively setting e± — and keep- 
ing /</0we obtain 



1 



1 



Ml I 2 



which exhibits a discontinuity as /i — > 0. The case F = 
was studied by Read and GreenSi who suggested that the 
fx = limit is identified with a critical point for which 
the ground state is the Haldane-Rezayi stated However 



7 



the situation is much the same for generic F. Let Pm 
denote the projection operator onto the subpsace for M 
Cooper pairs. If we take the projected mean-field ground 
state at /j, — (with e\ ^ 0) we obtain 

( \ M 
jyA |0) (38) 

for some normalisation constant C. Eq. (|38[) is exactly 
the dressing of the vacuum state that was discussed in 
Sect. HH and appeared in Sect, [ill] in the context of 
the Bethe ansatz solution. Thus the projected mean- 
field ground state (|55|) for the critical point at which the 
excitation spectrum is gapless (/i = 0) is seen to have a 
universal character in that it is independent of both F 
and G. 



V. CONCLUSION 

Using the conserved operators and Bethe ansatz equa- 
tions of the Richardson-Gaudin pairing Hamiltonian we 
obtained an integrable, exactly solvable, extended d + id- 
wave system. This model is curious in that the energy ex- 
pression depends quadratically on the roots of the Bethe 



ansatz equations, and also in that the Bethe ansatz equa- 
tions involve terms which do not allow for an electrostatic 
analogy. Nonetheless we were still able to study this sys- 
tem by solving the singular integral equation obtained 
from the Bethe ansatz equations in the continuum limit. 
We also showed that these results and those obtained 
in a mean-field approximation coincide. For a general 
Hamiltonian with two coupling parameters (which in- 
cludes the integrable case) it was found that the projected 
mean-field ground-state is universal at the critical point 
of gapless excitations. This suggests that results from 
the exactly solvable extended d + id-wave system may 
prove fruitful in understanding the d + id-wave system. 
A precedent for such an approach has been set in the 
different context of spin ladder systems, whereby it has 
been shown that extended exactly solvable Hamiltonians 
can provide extremely useful tools for making compar- 
isons between theoretical calculations and experimental 
data^l 
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